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ABSTRACT 


The  Kalman  Filter  is  implemented  to  process  range  data  output  from 
the  Cubic  Model  40  Autotape  system,  a  surface  position  locating  system 
currently  employed  on  the  underwater  tracking  ranges  at  Dabob  Bay  and 
Nanoose.  Results  are  presented  for  different  measurement  noise  and 
forcing  function  noise  statistics. 
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I.   INTRODUCTION 

The  Cubic  CM-40  Autotape  is  a  microwave  distance  measuring  system 
used  (by  the  U.S.  Navy  at  its  acoustic  underwater  tracking  ranges  at 
Dabob  Bay  and  Nanoose)  to  provide  reference  position  information  for 
units  on  the  surface  and  in  the  air  above  the  range.  This  portable 
system  consists  basically  of  an  interrogator  which  is  operated  aboard 
the  unit  to  be  tracked,  two  responders  operated  at  two  different  shore 
sites  and  the  associated  antenna/RF  assemblies.  Required  support  systems 
include  a  data  display  and  recording  setup  and  an  ADP  facility  for  off- 
line processing  of  the  Autotape  data.  Figure  1  shows  the  Autotape  system 
components  and  Figure  2  shows  a  typical  application  geometry. 

Historically,  the  Autotape  has  been  used  in  such  applications  as 
tracking  hydrophone  array  survey,  buoy  and  hydrophone  array  planting  and 
as  a  reference  position  indicator  for  calibrating  other  position-finding 
devices  against.  Generally,  the  Autotape  has  been  used  where  an  extremely 
high  degree  of  accuracy  is  not  required. 

In  operation,  the  system  will  provide  for  the  display  and  recording 
of  two  ranges  simultaneously,  once  per  second,  the  ranges  being  those 
between  the  interrogator  and  each  of  the  responders.  The  ranges  are 
computed  from  the  phase  delay  between  the  output  of  the  modulation  signal 
generator  and  a  signal  which  has  traveled  from  the  interrogator  to  a 
responder  and  back.  Ranging  accuracy  is  stated  by  the  manufacturer  to 
be  +  0.5  meter  +  10  ppm  x  range.  Ranging  frequencies  of  1500  KHZ,  150 
KHZ  and  165  KHZ  modulate  a  3000  MHZ  carrier,  yielding  a  maximum  unambiguous 
range  of  10,000  meters  with  a  resolution  of  0.1  meter.  However,  independent 
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FIGURE  1:     Cubic  Model   40  Autotape  System 
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testing  by  the  U.S.  Navy  [Reference  1]  has  shown  that  system  accuracy  may 
not  be  quite  as  good  as  stated  by  the  manufacturer. 

The  accuracy  of  the  Autotape  system  is  principally  dependent  upon 
range  errors,  the  geometry  of  the  system  and  the  method  of  data  reduction. 
These  factors  are,  in  turn,  affected  by  propagation  velocity,  system 
stability,  range  dependency,  land  survey  accuracy,  system  geometry, 
slope  reduction  and  data  smoothing.  A  final  anomaly  which,  depending 
upon  the  application,  can  substantially  degrade  the  quality  of  the  data- 
stream  out  is  the  orientation,  over  time,  of  the  interrogator  antenna  in 
the  vertical  dimension.  The  interrogator  antenna  has  only  a  10  degree 
vertical  beam  width.  Thus,  if  the  system  is  being  used  on  a  platform 
such  as  a  moderately  maneuvering  helicopter  or  a  ship  rolling  substan- 
tially in  the  seaway,  the  system  tends  to  frequently  lose  track,  resulting 
in  fairly  long  streams  of  useless  data. 

Present  data  reduction  techniques  employed  when  the  system  is  used 
on  either  of  the  ranges  (Dabob  or  Nanoose)  employ  two  overall  iterations. 
The  first,  or  initial  processing,  administers  the  following  three  correc- 
tions to  the  raw  range  data: 

1.  Range  Calibration  Correction:  This  is  a  fixed  value  (meters)  added 
to  or  subtracted  from  each  range. 

2.  Propagation  Velocity  Correction:  This  is  a  variable  correction  due 
to  the  atmospheric  index  of  refraction  at  the  particular  time  and 
place  of  the  exercise. 

3.  Slope  Reduction  Correction:  This  reduces  both  range  measurements 
(which  are  actually  slant  ranges  because  the  interrogator  and  the 
responders  are  not  normally  located  at  the  exact  same  elevation) 
to  a  common  horizontal  plane  at  sea  level. 

Subsequent  processing  of  the  data  includes  conversion  of  the  corrected 

ranges  to  a  rectangular  x-y  range  coordinate  system  and  a  moving  average 

smoothing  technique  which  employes  curve  fitting  algorithms  (linear, 


parabolic  or  logarithmic)  to  reduce  the  data  to  its  final  form.  Not 
uncommonly,  as  a  result  of  the  total  reduction  effort,  the  net  remainder 
is  an  inadequate  data  package  (in  terms  of  quantity)  for  proper  final 
evaluation. 

Figure  3  is  a  rectangular  plot  of  the  raw  ranges  recorded  during  a 
recent  array  survey.  The  purpose  of  this  project  has  been  to  design  a 
filter,  a  Kalman  filter,  which  would  provide  more  accurate  range  data, 
as  well  as  one  that  would  track  through  the  periods  of  "lost  track" 
ranging,  thereby  providing  a  significantly  larger  final  volume  of  data 
for  evaluation.  This  paper  presents  the  basic  theory  necessary  and 
includes  the  final  version  of  the  filter. 
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FIGURE  3:  Rectangular  Plot  of  Raw  Range  Data 
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II.  THE  FILTER  THEORY  AND  DESIGN 

A,  THE  SYSTEM  DYNAMIC  MODEL 

A  common  application  for  the  Autotape  system  is  its  use  as  a  reference 
position  locator  on  the  surface  unit  conducting  an  acoustic  hydrophone 
array  (range)  survey.  The  usual  exercise  plan  will  call  for  a  service 
unit,  carrying  the  interrogator  and  equipped  with  an  acoustic  pinger 
mounted  on  the  underwater  hull,  to  transit  three  concentric  circular 
tracks,  centered  above  the  array,  with  track  radii  ranging  from  100  to 
1,000  meters,  at  speeds  of  up  to  eight  knots.  The  direction  of  rotation 
for  the  outer  track  will  normally  be  opposite  to  that  of  the  middle 
circle.  While  the  service  unit  is  being  tracked  via  Autotape,  it  is  also 
being  tracked  by  the  acoustic  array.  By  comparing  the  acoustic  position 
data  with  that  from  the  Autotape,  a  digital  computer  is  able  to  compute 
actual  position  and  attitude  of  the  array. 

The  desired  estimates  will  be  those  of  position  and  velocity,  Rl,  R2, 

Rl,  R2.  It  is  proper  at'  this  point  to  define  a  number  of  terms  and  to 

summarize  some  pertinent  results  of  observer  theory.  First,  we  may  define 

a  fourth  order  state  vector: 

Rl 
R2 
Rl 
R2 

Recall  that  a  linear  system  can  be  described  in  the  continuous  time 

domain  as: 

x  (t)  =  A  x  (t)  +  D  w  (t) 


x  = 
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where:       x(t)   is  the  n-element  column  vector  of  the  states 

A  and  D^  are  nxn  and  nxp  matrices  describing  system  dynamics 

w(t)  is  a  q-element  vector  of  random  noise  inputs  to  the 
system 

The  system  measurements  may  be  expressed  as: 
z(t)  =  H  x  (t)  +  v  (t) 

where:        z_(t)  is  the  q-element  vector  of  system  measurements 
H.  is  the  qxn  weighting  matrix  for  the  measurements 
y_(t)  is  the  q-element  vector  of  random  measurement  noise 

The  corresponding  linear  discrete  model  may  be  written  as: 

x(k  +  1)  =  0  x(k)  +  r  w(k) 
with  no  deterministic  inputs  to  the  system. 


Also, 


z(k)  =  H  x(k)  +  v(k) 


For  the  system  under  consideration,  it  can  be  shown  that  the  state 
transition  matrix 
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for  a  sampling  interval  T  of  1  second.  A  block  diagram  of  the  system  is 
shown  in  Figure  4. 
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FIGURE  4:  Block  Diagram  of  Discrete  Linear  Estimator 
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The  following  assumptions  will  be  made  regarding  the  noise  processes 
and  the  initial  state,  x(o)  of  the  plant  [Ref.  2]: 

The  measurement  noise  has  zero  mean,  is  uncorrected,  and 

E  [  v(k)  v  (j)]  s  R(k)  ?..,  where  2  is  the  kronecker  delta 

The  forcing  noise  has  zero  mean,  is  uncorrelated,  and 
E  [  w(k)  wT(j)]  =  Q(k)  £kj 

The  forcing  noise  and  measurement  noise  are  uncorrelated. 

The  initial  state  is  a  random  variable  with  known  mean  and  covariance, 
and 

E[$x(o)-x0}   (x(o)  -  *0)   T]  =  PQ 

The  measurement  noise  and  initial  state  are  uncorrelated. 
The  forcing  noise  and  initial  state  are  uncorrelated. 

The  Kalman  Filter  equations  and  their  derivation  are  well  known 
[Ref.  2],  [Ref.  3]: 

G(k)  =  P(k/k-l)  HT(k)  [H(k)  P(k/k-l)  HT(k)  +  R(k)]"1  (1) 

P(k/k-l)  =  0  P(k-l/k-l)  0T  +  Q  r   (2) 

P(k/k)  =  [1  -  G(k)  H(k)]  P(k/k-l)  ^   (3) 

<  x(k/k)  =  x(k/k-l)  +  G(k)  [z(k)  -  H(k)  x(k/k-l)]  (4) 

x(k/k-l)  =  0(k/k-l)  x(k-l/k-l)  +  r(k/k-l)  w(k-l)  (5) 

* 

Where  the  notation  (k/k-1)  interprets  as  the  value  of  the  parameter  of 
note  at  time  k  given  measurements  at  times  up  to  and  including  time  k-1. 
(k/k)  and  (k-l/k-1)  have  similar  interpretations.  The  x_  denotes  the 
estimate  of  x.- 
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G(k)  represents  the  filter  gain  at  time  k.  P_  represents  the  covariance 
of  estimation  error; 


P  (k/k)  =  E  [e(k/k)  e'(k/k)]  =  E 


/h 


=  E  { 


ej(k/k) 


e^k/k) 
e2(k/k) 


en(k/k) 


[e^k/k)  e2(k/k)- 


-en(k/k)] 


e^k/k)  e2(k/k)  —  e^k/k)  en(k/k) 


e2(k/k)  e^k/k)   e2(k/k) 


—  e2(k/k)  en(k/k) 


\fc 


en(k/k)  e^k/k)   en(k/k)  e2(k/k)  —  e^(k/k) 


where  e_(k/k)  =  S(k/k)  -  _x(k).  A  complete  standard  block  diagram  for  the 
filter  and  an  information  flow  diagram  are  included  as  Figures  5  and  6 
as  slightly  different  viewpoints  from  which  the  system  may  be  viewed  and 
understood.  Figure  7  shows  a  timing  diagram  of  the  various  quantities 
contained  in  the  filter  equations. 
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FIGURE  5:  Kalman  Filter  Block  Diagram 


17 


H(k-l)        R(k-l) 


1 


Compute 
P(k/k)  =  [I-G(k)H(k)]  P(k/k-l) 


Q(k-l)      0(k-l) 

J ±_ 


Compute 
P(k/k-l)  =  0  P(k-l/k-l)0T+  Q 


i 


set  k 


H(k) 


R(k) 


Compute 
G(k)  =  P(k/k-l)HT  [H(k)P(k/k-l)HT(k)  +  R(k)]"1 


z(k) 


G(k) 


w(k-l)  0(k-l)  r(k-l) 


Compute 

x(k/k)  =  x(k/k-l)+G(k) 
[z(k)-H(k)x(k/k-l)] 


set  k 


*. ± 


Compute 


x(k/k-l)  =  0x(k-l/k-l) 


x(k/k-l) 


FIGURE  6:  Simplified  Information  Flow  Diagram  of  a  Discrete  Kalman  Filter 
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FIGURE  7:  Timing  Diagram  of  Filter  Equation  Quantities 
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B.  THE  PROCESSOR 

Appendix  A  is  a  flowchart  of  the  Kalman  filter  program  utilized. 
Initially,  the  matrices  describing  the  physical  system,  the  noise 
statistics  and  other  program  parameters  are  read  into  storage  and  printed 
out.  The  discrete  state-transition  matrix,  Phi,  is  computed  and  printed 
out  and  the  gain  schedule  is  computed  and  printed  out.  It  is  seen  that 
the  elements  of  the  gain  matrix  reach  a  steady  state,  and,  for  example, 
with  both  the  R  and  (^matrices  being  identity  matrices,  the  gain  reaches 
steady  state  between  k=5  and  k=10.  Therefore,  in  the  main  iteration  loop, 
the  filter  will  essentially  be  a  constant  gain  filter  for  k>  10. 

Next,  the  main  iteration  loop  commences.  The  initial  measurements 
are  read  and  xl(0/-l)  and  x2(0/-l)  are  initialized  to  these  values. 
x3(0/-l)  and  x4(0/-l),  representing  the  rates,  are  set  to  the  mean 
constant  value  (in  the  respective  directions)  of  4.0  meters  per  second. 
The  Autotape  output  is  a  5  significant  figure  output,  modulo  10,000, 
reading  to  0.1  meter.  Inherent  in  the  output  is  a  major  degree  of  jitter 
in  the  two  most  significant  digits,  which  would  significantly  distort 
the  covariance  of  measurement  noise.  Therefore,  as  an  option,  measure- 
ments could  be  gated,  and  the  gain  automatically  set  to  zero  in  those 
cases  where  the  residue  falls  outside  of  a  maximum  reasonable  bound. 

Commencing  with  k=0,  and  utilizing  the  known  values  for  £(0/-l)  and 
PJ0/-1),  the  Kalman  filter  equations  are  solved  iteratively  in  the 
following  manner  [see  page  15,  equations  (l)-(5)]: 

(1),  (3),  (4), 

Increment  k  to  k=l 

(5),  (2),  (1),  (3),  (4), 

Increment  k  to  k=2 

(5),  (2),  (1),  (3),  (4), 

etc. 
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Also  computed  on  each  iteration  are  the  error  residues: 

RES  =  z  -  x(k/k-l) 
and  the  one-step  prediction  errors: 

ERR  =  x(k/k)  -  x(k/k-l) 
Finally,  the  computations  are  tabulated  and  plots  are  produced. 

C.  NOISE  AND  ERROR  CONSIDERATIONS 

Reference  1  documents  an  Autotape  evaluation  which  was  conducted  in 
1971.  The  error  geometry  is  shown  in  Figure  8.  Graphically,  position 
is  determined  by  locating  the  crossing  point  of  the  two  range  arcs,  in 
conjunction  with  a  knowledge  of  the  baseline  formed  by  the  two  responders 
Since  each  range  has  an  associated  standard  deviation  (error),  the  point 
can  actually  be  enclosed  in  a  parallelogram  which  defines  the  probable 
position  within  one  standard  deviation  of  the  ranges.  The  shape  of  the 
parallelogram  will  vary  with  the  position  of  the  crossing  point  relative 
to  the  baseline,  as  indicated  in  Figure  8.   It  can  be  shown  that  the 
maximum  probable  error  (MPE)  will  be  minimized  where  the  range  arcs  are 
orthogonal.  Figure  9  diagrams  error  contours  which  are  actually  the 
locii  of  constant  MPE  for  two  particular  responder  sites  on  the  Nanoose 
Range.  Table  1  summarizes  pertinent  results  of  the  study. 
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TABLE  1 

Average  Range  Errors  (feet) 

Survey 

No.  Points 

R-l 

R- 

■2 

Error    Standard 
Average   Deviation 

Error 
Average 

Standard 
Deviation 

Array  04 

30 

-  0.5       2.8 

-  0.1 

2.8 

Array  07 

49 

-  1.3       2.3 

-  0.4 

2.4 

Array  08 

10 

-  0.8       4.4 

1.6 

2.8 

Array  09 

25 

3.8       2.6 

0. 

2.2 

Average 

0.3       3.0 

-  0.5 

2.6 
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For  the  purpose  of  modeling  the  covariance  of  excitation  noise,  it 
was  assumed  that  the  service  unit  transited  an  800  meter  circle  at  an 
average  speed  of  eight  knots.  Then: 

2 


a  = 


R 


(8  kts) 


[1830  meter] 
t  n.   mile    1 


3600 


sec 
Hr 


0207 


800  meters 
m 


sec 


Filter  performance  was  investigated  for  Q.  =  I_,  .11.,  and  .0H_, 


for 


and 


P(0/-1)  =  Pq  =  E 


a-  [!!] 


(fcW  -  *]  |i(0)  -  x]J 


7 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

L° 

0 

0 

1 

where  the  a  priori  x_(0/-l)  is  known  to  be  a  reasonably  good  estimate  -- 
approximately  the  same  accuracy  as  an  observation. 

D.   PROCESSOR  PERFORMANCE;  AUTHOR'S  CONCLUSIONS 

Table  2  summarizes  a  comparison  of  the  Kalman  filter  performance  with 
the  results  of  the  (corrected)  processing  by  the  program  presently  being 
used  for  the  cases  0.  =  1,  R  =  I>  and  Q  =  0.H,  R  =  1-  Figures  10,  11, 
12  and  13  are  residue  and  error  plots  for  the  example  Q_  =  .0H_,  R  -  I_. 

It  is  seen  that  the  Kalman  filter  will  satisfactorily  handle  the  data 
where  the  measurement  noise  statistics  approximate  those  used  in  the 
model.  However,  for  the  noise  resulting  from  the  jitter  which  appears 
in  the  "hundreds"  and  "thousands"  digits,  the  filter,  as  configured 
without  a  gate,  will  estimate  with  considerable  error.  The  raw  range 
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R2  was  clean  of  this  particular  noise  element,  and  the  results  as 
indicated  by  Figures  12  and  13  were  superior  to  those  for  Rl. 

It  is  suggested  that  the  Kalman  filter  be  used  as  the  first  iteration 
processing  of  the  Autotape  output. 
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FIGURE  10:  Residue  1  vs.  Time.   Q  =  .011,  R  =  I. 
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FIGURE  12:     Error  1  vs.   Time.       Q  =   .011,   R  =   I 


30 


o 
o 

CD 


O 
Q 

CD 


a 

Q 


o 
a 


CN 

i 


o 

a 


o 

CD  + 

I 


o 
o 

w 

CO 

1 


^^lf'      ^*-*' 


-^- — <x. 


■/— v 


3.00  1.00 

TIME  KCX101 


5.00 


FIGURE  13:     Error  2  vs.   Time.       Q  =   .011,   R  =   I 
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III.   FUTURE  FILTER  IMPROVEMENTS 

The  filter,  as  designed,  will  process  by  off-line  (forward)  filtering 
of  the  range  measurements.  It  is  suggested  that,  as  an  effort  to  further 
improve  upon  the  quality  of  the  processed  data,  a  fixed-interval  smoothing 
algorithm  (the  initial  and  final  times,  0  and  T,  are  fixed,  and  the 
estimate^  (t/T)  is  sought)  be  incorporated. 

For  the  system  and  measurements  described  by: 
*_  =  Fx  +  Gw 
z   =  Hx  +  y_ 

the  equations  defining  the  forward  filter  are,  in  the  time  domain  [Ref.3]: 

£  =  fx  +  PhV1  [z  -  Hx] ,  x  =  x^  (1) 

A  =  f EL  +  PLT  +  GQGT  -  PhV^P,   P_(0)  =  P^  (2) 

m  Y         (\  Y 

To  write  the  backward  filter  equations,  setT  =  T  -  t.     Then  —ms-  =  ~dt   , 


and 


Also, 


dx 

-p—  =  -Fx  -Gw,   for  0  <  T  <  T,  denoting  differentiation  with 

respect  to  backward  time. 


z(T)  =  Hx  +  v 


Then,  by  analogy,  the  backward  filter  equations  can  be  written  by 
changing  _F  to  -£  and  G  to  -G,  resulting  in: 

t    =  -fl  +  P.hV1  [z  -  Hx!  ] 


and      "4r~-b  =  ~--Pb  "  -PbF-T  +  -^-    "  ?bH--l-?b  (3) 
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FIGURE  14:  Relationship  of  Forward  and  Backward  Filters 


backward  filter    -a   „ 
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*.  p 


forward  filter 


From  Figure  14,  it  can  be  seen  that  the  smoothed  estimate  at  time=T 
must  be  the  same  as  the  forward  filter  estimate  at  that  point,  i.e., 

x  (T/T)  =  t   (T) 
and      P  (T/T)  =  P  (T) 

which  yields  the  required  boundary  condition  on  PT  , 

P^1  (t=T)  =  0,  or  P^1  (T=0)  =  0  (4) 

with  the  boundary  condition  on  £  (T)  not  yet  known.  Therefore,  define 
the  new  variable: 

s(t)  =  P^1  (t)  ^  (t)  (5) 

and  since  'x,    (T)  is  finite,  it  follows  that: 

s  (t=T)  =0,  or  s  (T=0)  =  0.  (6) 
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_1 

Reformulation  in  terms  of  P.  yields 


d      p-1.     .p-1        dp        p-1 


dT  M)  4)      |    dT  4      4) 

Thus,  equation  (3)  can  be  written  as: 


p"1  =     P"1   F  +  FT  PT1  -  PT1  G  Q  GT  PT1  +  HT  R_1  H  (7) 


for  which  equation   (4)   is  the  appropriate  boundary  condition. 

Differentiating  equation  (5)  with  respect  to  T,  and  with  some 
substitution  and  manipulation,  we.  arrive  at: 

H?rl=    (^  -P^GQG1  I  i  +  ^R"1!  (8) 

for  which  equation  (6)  is  the  appropriate  boundary  condition.  Equations 
(1),  (2),  (7)  and  (8),  along  with: 

P"1  (t/T)  =  P"1  (t)  +  P^1  (t) 

x  (t/T)  =  P  (t/T)  [r1  (t)  $  (t)  +  P^1  (t)  x^  (t)] 

define  the  optimal  smoother. 

Many  forms  of  the  smoothing  equations  may  be  derived.  The  form 
proposed  for  use  in  this  particular  case  is  the  Rauch-Tung-Striebel  form, 
with  the  discrete-time  expressions  summarized  as  follows: 
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Smoothed  State  Estimate   £(k/N)  =  £(k/k)  +  /^  [£(k+l/N)  -£(k+l/k)] 

where 

^  =  P(k/k)  0(k)T  P(k+l/k)_1 
for  k  =  N-l 

P(k/N)  =  P(k/k)  +  A.  [P(k+1/N)  -  P(k+l/k)]  A^ 
also  for  k  =  N-l 


Error  Covariance 
Matrix  Propagation 


Solution  of  the  equations  would  proceed  as  follows:  As  an  example 
and  because  it  is  slightly  easier  to  see  when  actual  times  are  used, 
suppose  NN  =  100.  On  the  forward  filter  pass,  the  values  of  x_(k/k), 
£(k/k-l),  P_(k/k)  and  P_(k/k-l)  would  be  computed  and  stored.  On  the 
final  iteration  of  the  forward  pass,  with  K  =  NN  =  100, 

£(100/100)  =£(100/99)  +G(100)  [z(100)  -  H  £(100/99)] 

i.e.,  we  have  computed  and  stored  £(100/100) . 

Now,  the  smoothing  process  commences  in  the  reverse  direction. 
Decrement  k  to  k  =  NN-1  =  99,  then 

£(99/100)  =£(99/99)  +  A(99)  [£(100/100)  -  x( 100/99)] 
stored  stored      stored 


and  A(99)  =  P(99/99)  0T  PU00/99)"1 

stored  stored 
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let  k  =  NN-2  =  98,  then 

£(98/100)  =  £(98/98)  +  A(98)  [£(99/100)  -  £( 99/98)] 

stored         computed    stored 

last 
iteration 

and      A (98)  =  P( 98/98)  0T  PJ99/98)"1 

stored     stored 

Also,  for  each  of  the  two  preceding  iterations, 

P(99/100)  =  P( 99/99)  +  A(99)  [P( 100/100)  -  P.(  100/99)]  AT(99) 
stored   computed  stored      sjtored 

P( 98/100)  =  P(98/98)  +  A(98)  [P( 99/100)  -  P( 99/98)]  AT  (98) 
stored   computed  computed    stored  computed 

etc. 

It  is  seen  that  the  smoothing  process  does  not  involve  the  processing 
of  actual  measurement  data.  It  does,  however,  utilize  the  complete 
filtering  solution,  and  so  fixed  interval  smoothing  cannot  be  done  real- 
time, on-line.  It  must  be  done  after  all  the  measurement  data  are 
collected.  Consequently,  computation  speed  will  not  be  the  most  important 
factor.  Storage  requirements  could,  however,  conceivably  be,  in  that  the 
quantities  to  be  stored  on  the  forward  pass  are  arrays.  It  is  seen  that, 
should  an  exercise  run  in  excess  of  30  minutes,  retention  of  the  data  at 
each  mark  could  require  in  excess  of  100K  bytes  of  memory,  which  could 
limit  the  facilities  upon  which  the  processor  could  be  utilized. 
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APPENDIX  A:  Processor  Flowchart  Main  Program 


C      START    J 

] 

r 

Read 

N,  M,  ND,  MD,  LD, 
NN,  DT 

V 

Wri  te 
N,  M,  ND,  MD,  LD, 
NN,  DT 

Read 
R,  Q,  P(0/-1),  A,  D 

i 

i 

Write 
R,  Q,  P(0/-1),  A,  D 

Call 
PHIDEL 

1 

t 

Wri 

te 
r 
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Do  10  k=0,  20 

1 

' 

Call 

GAIN 

t 

W 
k,  G(k) 

rite 

,  P(k/k-l) 

Do  15  k=0,  NN 

\ 

f 

IV3  =  6 

' 

> 

Call  IOW 

YES 


© 
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Decode 
Itime,  z(2,l),  z(l,l) 


Call  PROD 
x(k/k-l)  =  0x(k-l/k-l) 


Call  PROD 
zCOR  =  H  x(k/k-l) 


Call  SUB 
zMOD  =  z-zCOR 
=  z-H  x(k/k-l) 


RES1  =  zlsl(k)  -  xia(k/k-l) 
RES2  =  z2)1(k)  -  x2  1(k/k-l) 
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ARES1  =  /RES1/ 
ARES2  =  /RES2/ 


Call  GAIN 


Call  PROD 
xCOR  =  G  zMOD 
=  G(z-H  x(k/k-l) 


0 


© 


>0 


>o 


i_4 


x(k/k)  =  x(k/k-l) 


0 


0 
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Call  ADD 

x(k/k)  =  x(k/k-l)  +  xCOR 

x(k/k-l)  +  G(z-H  x(k/k-l)) 


ERR1 
ERR2 


Xjdc/k) 
x2(k/k) 


x^k/k-1) 
x2(k/k-l) 


Write  8  (Bulk  storage) 
I time,  RES1,  RES2,  ERR1 ,  ERR2 


PRINT 
zl,  z2,  xl(k/k),  x2(k/k), 
RES1,  RES2,  ERR1,  ERR2 


0         (£) 
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=  0 


106V 


>0 


105 


Call   TRANS 
0  — »0T 


Call  PROD 
TEMP3  =  P(k/k)0 


T 


Call  PROD 
TEMP4  =  &  TEMP3 
=  0  P(k/k)  0T 


Call  ADD 
P(k/k-l)  =  TJEMP4  +  Q 
=  0  P(k/k)  0T  +  Q 
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© 


Call  TRANS 
H  — >HT 


Call  PROD 
TEMP  =  P(k/k-l)H 


T 


Call  PROD 
TEMPI  =  H  TEMP 
=  H  P(k/k-l)HT 


Call  ADD 
TEMPI  =  TEMPI  +  R 
=  H  P(k/k-l)HT  +  R 


Call  RECIP 
TEMP2  =  TEMPI 


-1 


=  [h  P(k/k-l)HT  +  r\ 


-1 


Call  PROD 
G  =  TEMP  TEMP2 
P_(k/k-l)HT[7[  P(k/k-l)HT  +  p7 


-1 
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0 


Call  PROD 
TEMP3  =  G  H 


Call  ADD 
TEMP3  =  H]_  +  TEMP3 
=  I  -  G  H 


Call  PROD 
P(k/k)  =  TEMP3  P_(k/k-l) 
=  [1  "  £  3  P_( k/k-1) 


C 


RETURN 


3 
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COMPUTER  OUTPUT 
Q  =  I ,     R  =■  I 


K 

RAW    ' 

->1 

,3  AW     . 

?2 

c  tj    t  r  p  r  0 

C 

4222 

» 4 

42S2 

.  2 

4C22.4 

i 

4523 

,  ? 

4373 

,  ]_ 

4S27.4 

im 

4C33 

*-> 

4372 

7 

4  U  ~  2  •  1 

T 

4  3  3  5 

,  2 

4372 

.4 

4  535  .3 

4 

4C33 

T 

4  35  3 
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Q  =  I,  R  =  I 
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Q  =  o.n,    R  =  i 
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2=  o.on,  r  =  i 
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4656.6 
4660-8 
4657.7 
4666-5 
4673.6 
4674-3 
4676.8 
4601-0 
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1C5624 
1C5625 
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105528 
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